#include "base.h"

int main(int argc,char **argv)
{

//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@///

          //////////////////
          // Input   ///////
          /////////////////

	  int iOa2;
	  int imom;
          char dir[1024];



          ////// Parsing command line

          for(int i=1; i<argc; i++) {
            if(strncmp(argv[i],"-Nclust",7)==0)  sscanf(argv[i],"-Nclust=%d",&Nclust);
            if(strncmp(argv[i],"-dir",4)==0)  sscanf(argv[i],"-dir=%s",dir);
            if(strncmp(argv[i],"-T",2)==0)  sscanf(argv[i],"-T=%d",&T);
            if(strncmp(argv[i],"-L",2)==0)  sscanf(argv[i],"-L=%d",&L);
            if(strncmp(argv[i],"-iOa2",5)==0)  sscanf(argv[i],"-iOa2=%d",&iOa2);
	    if(strncmp(argv[i],"-imom",5)==0)  sscanf(argv[i],"-imom=%d",&imom);
          }

   /////////////////////////
  //// Momenta set
  //////////////////////////////
  char momentafile[1024]; sprintf(momentafile,"../INPUT/mom.dat");
  read_momenta_data(momentafile);
  Equivalent_Momenta(SinP2);
  double *nn=(double*)malloc(sizeof(double)*(Nmom));
  compute_nn(nn);

  // IO

   char file[1024];
   char Oa2print[1024];
   char dir_qav[1024];
   sprintf(dir_qav,"%s/QAV/",dir);
   whichOa2(iOa2,Oa2print);

    //////////////////////////////////
    // Read Zq 
     /////////////////////////////

   double *Zq=malloc(sizeof(double)*(Nclust+1));
   O55 *DpPC=malloc(sizeof(O55)*(Nclust+1));
   O55 *DpPC_inv=malloc(sizeof(O55)*(Nclust+1));
   O55 *ZpPC=malloc(sizeof(O55)*(Nclust+1));
//   O55 *DeltapPC=malloc(sizeof(O55)*(Nclust+1));


   double *aux_clust=malloc(sizeof(double)*(Nclust+1));


  sprintf(file,"Zq_%s",Oa2print);
  read_input_onemomenta(dir_qav,file,Zq,imom);

  for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++)
  {
		sprintf(file,"DpPC%d%d_%s",iop1+1,iop2+1,Oa2print);
		read_input_onemomenta(dir_qav,file,aux_clust,imom);
		for (int iclust=0; iclust<Nclust+1; iclust++) DpPC[iclust][iop1][iop2]=aux_clust[iclust];
  }



  // temporal for the inversion
//   double **DeltaT=malloc(sizeof(double*)*(5));
//   for (int iop1=0; iop1<5; iop1++) DeltaT[iop1]=malloc(sizeof(double)*(5));
   double **D=malloc(sizeof(double*)*(5));
   for (int iop1=0; iop1<5; iop1++) D[iop1]=malloc(sizeof(double)*(5));
   double **DT=malloc(sizeof(double*)*(5));
   for (int iop1=0; iop1<5; iop1++) DT[iop1]=malloc(sizeof(double)*(5));
   double **Dinv=malloc(sizeof(double*)*(5));
   for (int iop1=0; iop1<5; iop1++) Dinv[iop1]=malloc(sizeof(double)*(5));

   int iflag;


   for (int iclust=0; iclust<Nclust+1; iclust++)
   {
	///////////////////////////////////
	// With Delta
	/////////////////////////////////
	
	/*
	// 1. Compute Delta
	compute_delta( DpPC[iclust], DeltapPC[iclust]);
	// 2. Correct Dynamical matrix
	// D_corrected=D[I+Delta^T]
	for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++) DeltaT[iop1][iop2]=DeltapPC[iclust][iop2][iop1];
	for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++) {
		D[iop1][iop2]=DpPC[iclust][iop1][iop2];
		for (int iop=0; iop<5; iop++)	D[iop1][iop2]=D[iop1][iop2]+DpPC[iclust][iop1][iop]*DeltaT[iop][iop2];	
	}
	// 3. transpose
	for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++) DT[iop1][iop2]=D[iop2][iop1];
	// 4. invert
	inverse_real(DT,5,Dinv,iflag);
	if (iflag==1) {break; printf("D matrix is singular\n");}
	for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++) DpPC_inv[iclust][iop1][iop2]=Dinv[iop1][iop2];
	// 5. Z=Zq^2 * D^T^-1
	for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++) ZpPC[iclust][iop1][iop2]=Zq[iclust]*Zq[iclust]*DpPC_inv[iclust][iop1][iop2];

	*/
	//////////////////////////////////////////////////////
	// Why compute delta? I can invert directly the all matrix 
	// It must be diagonalized in blocks in the chiral limit
	// Away from the chiral limit the termos out of the diagonal are related with the old deltas (but they're not directly them)
	//////////////////////////////////////////////////////////
	// 1. D matrix
	for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++) D[iop1][iop2]=DpPC[iclust][iop1][iop2];
	// 2. transpose
	for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++) DT[iop1][iop2]=D[iop2][iop1];
	// 3. invert
        inverse_real(DT,5,Dinv,iflag);
        if (iflag==1) {break; printf("D matrix is singular\n");}
        for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++) DpPC_inv[iclust][iop1][iop2]=Dinv[iop1][iop2];
        // 5. Z=Zq^2 * D^T^-1
        for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++) ZpPC[iclust][iop1][iop2]=Zq[iclust]*Zq[iclust]*DpPC_inv[iclust][iop1][iop2];

	
   }


     ///////////////////////
     /// WRITING OUTPUT
     ///////////////////////

//     printf(" WRITING OUTPUT IN %s\n",dir);

   for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++)
   {
         sprintf(file,"ZpPC%d%d_%s",iop1+1,iop2+1,Oa2print);
         for (int iclust=0; iclust<Nclust+1; iclust++) aux_clust[iclust]=ZpPC[iclust][iop1][iop2];
	 write_output_onemomenta_real(dir_qav,file,aux_clust,SinP2_eq, imom);
   }

/*   for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++)
   {
         sprintf(file,"DeltapPC%d%d_%s",iop1+1,iop2+1,Oa2print);
         for (int iclust=0; iclust<Nclust+1; iclust++) aux_clust[iclust]=DeltapPC[iclust][iop1][iop2];
         write_output_onemomenta_real(dir_qav,file,aux_clust,SinP2_eq, imom);
   }*/

    free(Zq); free(DpPC); free(ZpPC); free(DpPC_inv);
    return 0;

}

